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ABSTRACT 


A general  one  dimensional  transient  heat  transfer  numerical  program  has 
been  developed  for  composite  building  constructions  with  arbitrary 
air  gap  locations.  The  complete  Fortran  language  program  as  used  on  the 
NBS  Unlvac  1108  computer  is  given.  A discussion  of  the  program  and  in- 
struction for  its  use  are  facilitated  by  the  aid  of  examples.  Numerical 
solutions  using  the  present  program  compare  favorably  with  experimental 
data  in  standard  fire  endurance  tests. 
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1 . 0 INTRODUCTION 


The  fire  performance  of  building  constructions  is  generally  evaluated 
by  a large  scale  laboratory  fire  endurance  test'  (ASTM  E 119)  in  which 
one  surface  is  exposed  to  fire,  controlled  according  to  a prescribed  in- 
creasing temperature  history  simulating  the  burnout  of  combustibles. 

Hie  fire  endurance  rating  of  the  construction  is  the  time  period  during 
which  it  withstands  the  fire  exposure  without  (a)  structural  failure, 

(b)  the  development  of  cracks  through  which  flames  can  cross,  or  (c)  the 
temperature  rise  on  the  unexposed  surface  exceeding  a prescribed  limit 
(250°F  rise  average,  325°F  rise  at  a single  point).  Where  the  failure 
criterion  is  due  to  heat  transmission  without  complications  due  to 
structural  or  physical  effects,  heat  transfer  analysis  should  provide  a 
means  for  prediction  and  design. 

A particular  aspect  of  the  fire  endurance  test  which  is  not  well  defined 
but  which  probably  plays  a significant  role  in  fire  performance,  is  the 
effect  of  mass  flow  (air  and  combustion  gases)  due  to  pressure  differences 
since  typical  building  constructions  consist  of  a series  of  composite 
layers  and  intermediate  air  layers,  a transient  heat  and  air  infiltration 
model  was  formulated.  Hie  program  is  particularly  suitable  for  evaluat- 
ing the  thermal  fire  endurance  of  building  constructions  where  various 
combinations  of  solid-to-solid  and  solld-to-air  contacts  are  encountered. 
For  each  solid  layer,  the  present  program  has  provisions  for  phase 
changes,  heat  generation  and  absorption,  and  thermal  property  variations 
commonly  found  in  building  materials.  Hirough  the  air  spaces  the  modes 
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of  heat  transfer  include  radiation  and  convection  with  temperature -depen 
dent  heat  transfer  coefficients  and  air  properties. 

A number  of  analog  and  numerical  programs  for  fire  endurance  evaluations 
had  been  in  existence  for  some  time  [1,  2,  3,  4j.  A more  flexible  and 
general  finite  difference  program  was  developed  by  Krokosky  as  recently 
as  1970  [,5J  • For  a review  of  the  existing  thermal  analysers  for  fire 
endurance  evaluation  one  is  thus  referred  to  and  5J. 

The  present  numerical  program  was  developed  to  Incorporate  into  fire 
endurance  analyses  the  following  features  which  are  desirable  and  yet 
not  readily  available  in  existing  programs: 

1.  Options  to  handle  heat  exposure  on  one  or  both  sides  of  structure 

2.  Heat  balances  in  air  spaces  to  allow  for  air  infiltration,  and 
heat  generation  and  absorption  in  air  spaces . 

3.  For  ease  of  application  to  building  structures  an  input  card, say 
lOllOlj is  sufficient  to  Instruct  the  computer  of  the  specified 
number  of  solid  layers  and  air  spaces  in  a given  problem. 

4.  Temperature -dependent  properties  and  known  chemical  heat  exchange 
of  various  building  materials  are  stored  in  a subroutine  and 
called  by  an  input  card, say  2331,  where  the  numbers  indicate 
coded  materials  stored  in  the  subroutine. 

5.  Duration  as  well  as  amount  of  known  chemical  heat  exchanges  can 
be  varied  in  any  material. 

2.0  GOVERNING  EQUATION  AND  BOUNDARY  CONDITIONS 

Hie  governing  equation  for  one  dimensional  transient  heat  flow  is  the 
well  known  heat  diffusion  equation.  Including  a term  for  internal  heat 
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generation,  this  can  be  expressed  as, 


— (k— ) + g = Pc' 


dx  Bk 


(1) 


where , 

T=  absolute  temperature  in  solid. 

X = coordinate  in  direction  of  heat  flow, 
k = heat  conduction  coefficient. 

g = time  rate  of  heat  generation  per  unit  volume  in  solid. 

P=  density  of  solid, 
c = specific  heat  capacity  of  solid. 

P=  time 

Boundary  Conditions  [^6,  7^  t 

Hie  following  types  of  boundary  conditions  need  to  be  evaluated  for  the 
general  problem; 


(a)  solid  to  air 

(b)  energy  balance  in  air  space 

(c)  air  to  solid 

(d)  solid  to  surrounding 

(e)  solid  to  symmetry  plane 

(f)  solid  to  solid 

(g)  furnace  gases  to  first  surface  la3rer. 


Figure  1 
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(a)  Solid  to  Air  (See  Figure  1) 


p^c^Ax^ 


BT. 


4, 


Ax. 


2 se  ’^1  Bx  ■ ^ ■''^1  2 

(b)  Energy  Balance  in  Air  Space  (See  Figure  1) 

"a  %a  Pa  |t  = ' ^2^ 

(3) 


(c)  Air  to  Solid  (See  Figure  1) 


Be 


^"^2  4 4 ^^2 

= 1^2  ^ + aei2  (TA  - T2)  +h2(TA  - T2)  +g2 


(4) 


In  the  above  equations  the  subscripts  1 and  2 indicate  solid  1 and 

1 


solid  2 respectively  as  shown  in  Figure  1,  and  where  ei2  = 


and  £x  j Sg  are  emissivity  of  surface  1 and  surface  2, 
T25  T^:  surface'  temperatures  as  shown  in  Figure  1 

Ax^,  AX2:  incremental  length  in  x direction 

ct:  Stefan  Boltzman  constant 


J-  -1 


Si  S3 


TA:  temperature  of  air 


h^  = K 


T - TA \ — 

1 1+  : convection  heat  transfer  coefficient  from  solid  to  air 


L:  characteristic  dimension  of  panel 

k^,  k2t  heat  conduction  coefficient  in  solid  1,  and  2 
®1’  ^2*  rate  of  heat  generation  or  absorption  per  unit  volume  in 

solid  1 and  2 
£,  : air  gap  spacing 

3 

c : specific  heat  capacity  of  air 

pa 


p : density  of  air 

(TA  - T \~ 

12  4 : convection  heat  transfer  coefficient  from  air  to  solid 

L / 
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m:  rate  of  mass  transfer  per  unit  area  due  to  pressure  difference 


g : rate  of  heat  generation  per  unit  volume  in  air  space  due  to 

a 

combustibles, 

K:  an  empirical  convection  heat  transfer  constant,  K=,27  for  vertical 

surfaces,  and  K=,38  for  horizontal  surfaces. 

Figure  2 - y Radiation 

t 

/ ^ 

- Solid  > T 

* — — — ^ 

t 


Convection 


(d)  Solid  to  Ambient  (See  Figure  2) 

pcAx  Bt  Bt  4 T ^ 4-c 

where  the  additional  variables  are. 


Ax 

2 


(6) 


T : ambient  temperature 

o , 

It  - T \-- 

h = Kl  olA-  : convection  heat  transfer  coefficient 


e:  emissivity  of  surface 

T:  surface  temperature  of  solid 

Figure  3 - 


/ 

/ 


/ 

/ 

Solid  ^ T 

/I  Air  Space 


-| 

T ^ 

/. 

^ Solid 


• ■ Symmetry  plane 

(e)  Solid  to  Symmetry  Plane  (See  Figure  3) 

By  symmetry,  temperature  on  both  solid  surfaces  that  face  each  other 

are  equal,  however  T - T , So  there  will  be  heat  transfer  from  solid 
’ s 

to  air.  At  the  interface  we  have. 


(8a) 
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At  the  air  space  we  have, 


i c P 
s pa  a 


3T 

s 

w~ 


- h (T  - 
s 


+ cre(T 


4 


(8b) 


where, 

T:  surface  temperature  of  solids 

T^:  temperature  of  air  in  symmetry  plane 

distance  from  solid  to  symmetry  plane 

convection  heat  transfer  coefficient  from  solid  to 
symmetry  plane. 


Figure  4 - 
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Solid  1 Solid  2 

(f ) Solid  to  Solid  (See  Figure  4) 

Consider  interface  between  Solid  1 and  Solid  2 as  shown  in  Figure  4, 
Let  and  be  the  interface  temperature  in  each  solid  respectively. 
A heat  balance  for  the  interface  node  can  be  expressed  as,  = T^, 


Pi  c, 


Axi 


+ PoC 


Ax2 


Ax 


1 


Ax. 


Bt 


-k. 


1 


+ k. 


'l"l  2 se  ' ^2"2  2 Be  ^1  2 ^ ^2  2 "l  Bx  ' "2  Bx 

where  subscript  indicates  conditions  in  solid  layer  1 or  2 . 


(5) 


Figure  5 


Radiation 


Convection 


Solid 


(g)  Furnace  Gases  to  First  Solid  Surface  (See  Figure  5) 

Consider  heat  transfer  from  furnace  gases  to  first  solid  surface.  The 
main  modes  of  heat  transfer  will  be  radiation  and  convection. 
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A good  approximation  of  ASTM  E 119  fire  curve  is  given  "by  the  following 
three  formulas,  (where  T = temperature  in  degrees  C,  and  0 = time  in 
minutes ) ; 

T = 9^0  °C  + 20°C  0 < 0 < 50  min 

(error  ± 2.%) 

T = 926°C  + O.T0°C  - 0.0131(120  - O)^  °C 

50  min  < 0 < 115  min 
(no  appreciable , error) 

11  ^in  < 0 < U80  min 
(linear  exact) 
first  solid  surface, 

- tJ)  + hp(Tp  - T^)  + g ^ (9) 

temperature , and 


FINITE  DIFFERENCE  EQUATIONS  ^ 

Consider  general  one-dimensional  heat  transfer  problem  containing  N 
solid  layers  and  m air  spaces  in  any  order.  To  facilitate  formulation 
and  discussion  let's  introduce  i as  the  running  index  for  air  spaces. 
In  Figure  6,  a general  multi-layer  configuration  is  shown,  where 
numbers  In  circles  indicate  "the  applicable  equation  at  the  given 
node  as  discussed  in  the  previous  section.  ” . . 


T = 926°C  + O.70°C 


Heat  balance  from  furnace  to 

im 

2 iO  ^ ^x  ^^F 


where  T^  is  absolute  furnace 

V ■ ^ 


3.0  PROBLEM  FORMULATION  AND 
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Solid 

1=1 


Radiation 


D O- 


Solid 

i 

■ O'  C 


Air 


© 0 © 


Solid 

i+1 


Solid 


Solid 

1=N 


T. 


'll 


'i+1 


Nn 


Convection 


X 


Figure  6.  A Composite  Wall  with  N Layers 
Note:  Number  in  circle  indicates  applicable  equation  as  discussed  in 

Section  2.0. 

To  solve  our  problem  numerically  with  the  nonlinear  heat  diffusion 
equation  and  associated  complex  system  of  boundary  conditions  we  shall 
use  forward  time  differencing  and  central  space  differencing  scheme. 
Furthermore,  we  must  use  three  subscripts. 

Let : 


T.  Temperature  of  ith  node  in  ith  solid 

i:  1...N  (number  of  solid  layers) 

j:  l...n  (number  of  nodes  on  solid  layer  ) 

k:  1 . . .m  (number  of  air  spaces) 

See  Figure  7 for  finite  differencing  network. 


Figure  7.  Sketch  of  Finite  Differencing  Network  in  a Solid  Layer 


2 ■ 
1 


Ax 
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Ax 

3 


Ax 


Ax  ! Ax 


ith  solid  layer 


1 ^ 

! 2 

i n 

Lj 


nth  node 
j=l . . .n 
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Furthermore , let 


TAj^:  Temperature  of  kth  air  spacing 

0 : Time 

AO:  Time  increment 


T 


i ,J 


Temperature  at  (O  + AO)  of  jth  node  in  ith  solid 


Si.:  Thickness  of  ith  solid 

1 

SiA,  : Thickness  of  kth  air  spacing 

K 

i. 

Ax.  = — ^ 

1 n - 1 


Difference  spacing  for  ith  solid 


g.  .:  rate  of  heat  generation  per  unit  volume 

k.  heat  conduction  coefficient 

O' ^ : heat  diffusion  coefficient 

ga  : rate  of  heat  generation  per  unit  volume  in  air  space 

K 

o 

m^:  rate  of  mass  transfer  per  \init  area 

Applying  our  finite  differencing  scheme  we  have  the  following  general 
finite  difference  expressions  corresponding  to  the  previously  discussed 
governing  equation  and  boundary  conditions;  [6,  7,  and  8] 


Governing  equation  ^ 

1 2 ®i i^^i 

hj  = i:  (h(j  - 1) " h(j . i)> " " eV 

1 1 1 


(1) 


, ..  ( Ax.  ) 

where  M.  = i 
1 


a . A0 
1 


From  equation  (l)  we  require  > 2 for  stability. 
Solid  to  Air, 


T.'  = T.  + §-(T.  , - T.  ) - R.[(T.  - (T._^,  - 

i»n  i,n  M.  i,n-l  1 ,n  i i ,n  i+l,l 

^ 5_  2 

H.  (T.  - TA,  + ®i,n  1 

k.  M. 

1 ,n  1 


(2: 
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where 


ZABcre 

= 7 

1 Ax . p . c . 
11  1 


c^:  specific  heat  of  ith  layer, 


'12 


: where  and  are  emmisivity  of  two  parallel 


solid  layers. 


H. 


2A9K 


i 0 .c ,Ax. (L)7 
^i  1 1 


L:  characteristic  dimension  of  panel. 

Air  Space  Heat  Balance 

’’A;  ' ^ '’k<A.n  - \<h,n  - 


+ 


where 


(3) 


A9 


-CA,  c p 
K pa^a 


P^:  air  density  function  of  temperature. 

c : specific  heat  of  air. 

pa 


q = ^9 

k - 


A9 


iA.  c p 
k pa^a 


c 0 
pa  a 


Air  to  Solid 


t: 


= T. 


+ R.r  (T.  ) - (T. 


i+1,1  "i+1,1  ■ ^^i-'"in'  ""i+1,1^  ^ ‘ M. 


i+1 


(T  - T ) 

1+1,1  i+1, 2^ 


+ H 


i+1 


(TA^ 


T. 


i+1,1 


5, 


(4) 
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R.,  M.  , H,  ,,  as  defined  before. 
1 i+1  1+1 

Solid  to  Ambient 


rp  ; 

N,n 


2 

= T + — (T 

N,n  M^^'N,n-l 


- T. 


N,n 


- T_  ) - 


N N,n 


- To> 


^N,n^ 


(6) 


where 


ZABce 


N 


n'^n 


emisivity  of  solid  (last  node  to  ambient). 

and  M„  as  defined  before 
N TSf,n  N 

T^:  ambient  temperature 

Solid  to  Symmetry  Plane 

2 n^^^N^ 

Tm  = T..  + ^(T.,  T - - T )4  + -r^—r: 

N,n  N,n  N,n“l  N,n  N N,n  s 

4 4 

- R-,(T._  - T ) 

N N,n  s 


(8a) 


T'  = T + P (T^,  _ T )l  + fT  - T ^) 

s s s N,n  s N N,n  s 


(8b) 


where 


: temperature  at  symmetry  plane 


P =1  K \ 4S_ 

= ' L^l 


I : distance  between  solid  surface  and  S5mimetry  plane, 

s 


Solid  to  Solid 


rp  ( _ rn  5 

i+1,1  i,n 


12 


= [A.T.  + A.  ,^T.  T + B.g,  + B. 

i,n  ^ 1 i,n  i+1  1+1,1  i 


’+l®i+l,l 


D.  (T.  - T.  ) - D.,,  ,(T.,,  , - „)](A.  + A_,  ) 

i,n  i,n-l  i,n  i+l,l  i+l,l  i+l,2  i i+l 

(5) 


where 


p.c .Ax. 

A = ^ ^ ^ 

i 2 


Ax.  A0 


k.  A9 


i , n Ax , 


i+1,1  Ax.,t 
1+1 


rate  of  heat  generation  per  unit  volume  in  ith  solid. 
Furnace  Gases  to  First  Solid  Surface 


' A,i  - A, 


^1  l^^l 

1 > ^«i^^F  - ^i,i> 


— (T  - T ) 

1,1  1,2^ 


(9) 


where 


2A9o'e^ 

A ' Ax,p^Cj^ 

e^t  emisivity  of  first  surface  (furnace  gases  to  1st  solid). 
TpC  furnace  temperature 

H = ^ 2A9 

1 p^c^Ax^ 
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4.0  DISCUSSION  OF  NUMERICAL  PROGRAM 


A flow  chart  of  the  main  program  is  presented  in  Figure  7.  A listing  of 
the  complete  numerical  program  is  also  presented  at  the  end  of  this 
report.  The  logical  sequence  of  the  main  program  can  be  grasped  read- 
ily by  first  considering  the  flow  chart. 

A list  of  input  parameters  is  read  in  on  data  cards  (examples  of  which 
are  shown  at  Section  5.3^as  follows: 

NN=  number  of  solid  layers 
N=  number  of  nodes  in  each  layer 
M=  number  of  air  spaces 

ID=  sequence  of  numbers  1 or  0,  specifying  the  sequence  of  solid  layers 
and  air  spaces,  e.g.,  101011101  means  solid-air-solid-air-solid- 
sol  id  -sol  id  -air  -sol  id  . 

IDD=  sequence  of  numbers  specifying  the  material  of  each  solid  layer, 
e.g.,  122321  means  the  six  solid  layers  of  the  problem  are  of 
materials  1,  2,  2,  3,  2,  1 in  that  order.  Various  temperature 
dependent  material  properties  are  stored  in  Subroutine  Prop. 

AL(I):  thickness  of  ith  solid  layer  in  feet 

ALPHA(I):  heat  diffusion  coefficient  of  ith  solid  layer,  ft^/hr. 

(ALPHA  assumed  constant  other  temperature  dependent  thermal  properties 
stored  in  Subroutine  Prop. 

RHO(I):  density  of  ith  solid  layer  in  Ib/ft^ 

AI(I):  thickness  of  air  gap  spacing  in  feet 

GA(k):  rate  of  heat  generation  per  unit  volume  due  to  combustibles  in 

kth  air  space,  Btu/ft^-hr. 

DTHETA:  time  increment  in  fraction  of  hour. 


H 


AMO:  Mass  flux  through  walls  in  Ib/ft^-hr 

GMl,  GM2,  GM3  etc:  Heat  absorption  or  release  due  to  phase  change  in 

material  1,  2,  or  3 etc.,  Btu/ft^hr. 

Subroutine  Prop:  stores  temperature  dependent  thermal  properties  and 

phase  change  reactions  of  some  common  building  materials.  This  sub- 
routine can  be  expanded  as  desire  when  new  materials  are  encountered. 
When  called  from  the  main  program  this  subroutine  supplies  the  thermal 
properties  for  the  ith  layer  of  solid  currently  being  calculated.  The 
function  subroutines  are  self  explanatory.  The  comment  statement  at  the 
beginning  of  each  function  subroutine  properly  identifies  it  with  the 
corresponding  equation  and  boundary  conditions  in  section  3.0. 
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5.0  SAMPLE  CASES  CALCULATED 


5.1  Description  of  Panels 

Case  1:  Four  solid  layers  (gypsum,  wood,  wood,  gypsum), 

heat  from  one  side  as  shown: 


Air 
3 


p5/8"  gypsum,  board 
Air 


Air 

^3"- 


-s-3"^ 


1/2”  plywood 


To 


three  air  gaps. 


Case  2:  Four  solid  layers  (gypsum,  wood,  wood,  gypsum),  three  air  gaps, 

heating  from  both  sides  as  shown: 


Case  3:  One  solid  layer  (plywood  either  1/2”  or  5/8”  thick),  heating 

from  one  side  of  panel. 


Case  4:  Two  solid  layers  and  one  air  gap  (Brick  - Air  - Brick), 

heating  from  one  side  as  shown: 

i ( 

I Brick  j Air  | Brick  I 

T,^  U-3  3/4'^~tf — 6”  5»K-3  3/4”-^  To 

" 1 ! ‘ ■ 

i ! 

5.2  Thermal  Properties 

The  following  thermal  properties  taken  from  Ref.  1 and  Ref.  2 are  used 
in  all  calculations,  e=,9  for  all  surfaces. 

Gypsum  Board 
O'  = .008  ft^/hr 
p = 60  Ib/ft^ 
c = .26  Btu/lb°F 
K = .125  Btu/hr  ft°F  0<T<200°F 
K = .075  Btu/hr  ft°F  200°F<T<400°F 
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K = (o05-t-  T ) Btu/hr  f t °F  400 °F<T<2000 °F 
16,000 

3 

Heat  of  moisture  desorption  and  calcination:  20,740  Btu/ft 


Plywood 

,006  ft^/hr 
p = 31  Ib/ft^ 

c = ,67  Btu/lb°F 
K = ,065  Btu/hr  £t°F 


Heat  of  moisture  desorption: 


15  00  Btu/ft^ 


Brick 

O' = ,028  ft^/hr 
p = 110  Ib/ft^ 

c = .216  Btu/lb°F 
K =1,0  Btu/hr  ft°F  0<T<200°F 

K = .46  + 2T/10,000  Btu/hr  ft°F  200°F<T<2000°F 

3 

Heat  of  moisture  desorption  = 5800  Btu/ft 
5.3  Input  Data 

The  following  are  print  outs  of  input  data  for  the  sample  cases 
calculated: 


Case  1: 

F our 

solid  layers. 

3 air 

gaps,  furnace  on 

one  side 

of  panel. 

(gypsum 

- air 

- wood  - air 

- wood 

- air  - gypsum) 

4 

5 

3 1 

68,000 

,000  5380,000 

1500.000 

1 0 1 

0 1 

0 1 

1 3 3 

1 

,250 

,000 

,250 

.000 

,250 

,000 

,052 

,008 

60,000 

,260 

,042 

,006 

31,000 

.670 

,042 

,006 

31,000 

,670 

,052 

,008 

60,000 

,260 

1 7 


Case  2: 

Four  solid 

layers,  3 air 

gaps  furnace  on 

both  sides  of  panel. 

(Note: 

Symmetry  ha 

s been  invoked 

gypsum  - air  - 

wood  - symmetry  plane) 

2 

5 1 

2 68.000 

.000 

5380.000  225O0OOO 

10  1 
1 3 

„500 

„052 

oOOO 

o008 

6O0OOO 

.260 

.042 

.006 

31.000 

.670 

Case  3: 

One  solid 

layer  (plywood) 

furnace  on  one 

side  of  panel. 

1 

1 

5 0 

1 68.000 

.00  5800.000 

1500.000 

3 

.042 

.006  31 

.000  c670 

Case  4: 

Two  solid 

layers  and  one 

air  gap  (brick 

- air  - brick)  furnace 

on  one  side  of  wall„ 


2 

2 5 1 

1 68.000 

.00 

5800.000  1500.000 

10  1 
2 2 

.5 

.000 

.313 

0O23 

110.000 

.216 

o313 

.023 

110.000 

.216 

Compari 

son  plots  of 

calculation  and 

standard 

fire  endurance  tests  are 

presented  in  figo  8,  9,  and  lOo  Tests  were  conducted  at  NBS,  Washington, 
D,C,  and  National  Gypsum  Corporation,  Buffalo,  N.Y. 


Finally  a word  of  caution  in  the  use  of  a thermal  analyzer  program  with 
air  gap  heat  balances  such  as  the  present  program.  Due  to  the  low  heat 
capacity  of  air,  one  must  exercise  caution  in  selecting  the  incremental 
A 0 andAXsuch  that  the  air  gap  temperature  will  not  rise  above  the 
temperature  of  the  node  preceeding  it.  From  experience  when  this  occurs 
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instability  will  set  in.  The  solution  is  then  to  decrease  A0  or  X 
and  at  the  same  time  maintaining  the  criteria  The  time  increment 

A0  is  5 seconds  and  the  number  of  nodes  is  5 in  all  calculated  cases 
except  case  4 where^dis  decreased  to  2,5  seconds  and  n is  increased  to 
25  to  avoid  the  instability  problem  discussed  in  above. 
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Fig.  8 Temperature  History  of  Case  1 (4  Layers, 3 Air  Gaps,  Gypsum-Wood-Wood-Gypsum,  Heating  from  one  side). 


do'3b'niVy3d(;J31 
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Fig.  9 Temperature  History  of  Case  2 (4  Layers,  3 Air  Gaps,  Gypsum-Wood-Wood-Gypsum, 

Heating  from  both  sides). 


V;  TF 
A T(UD 

O T(l,r-!)  (CALCULATLiD  r-Ofi  1/2"  PLYWOOD) 
0 T(i,n)  (F-X'PERiMbMTFOR  '/£"  PLYWOOD) 
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FIG  8 FLOV/  chart  of  TRANSIENT  HEAT  TRANSFER  THROUGH  COMPOSITE  WALLS  WITH  ARBITRARY 
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I jPUTS 

NN=  NUMBER  UF  SOlIDS 

N=  number  of  nodes  in  solid 

M = NUMBER  OF  AIR  SPACES 

ITE5T  = AN  integer  THAT  INDICATES  WHETHER  FORMULA  6 OR  FORMULA  8 
WILL  BE  USED  IN  COMPUTATIONS?  1 INDICATES  FORMULA  6 ANY  OTHER 
NUMBER  INDICATES  FORMULA  8 
TO  = INITIAL  TEMPERATURE  NOW  6P 

GMl  = HEAT  RELEASE  ^ER  UNIT  VOLUME  DUE  TO  MOISTURE  VAPORIZATION  FOR 
GYPSUM 

bM2  = HEAT  RELEASE  PER  UNIT  VOLUME  DUE  TO  MOISTURE  VAPORIZATION  FoR 
BRICK 

GM3  = heat  release  PER  UNIT  VOLUME  DUE  TO  MOISTURE  VAPORIZATION  FOR 
WOOD 

ID  AN  ARRAY  OF  INTEGERS  THAT  INDICATES  THE  SOLID-AIR  CONFIGURATION 
C THATS  BEING  TESTED  0=  AN  A I R GAP  NON-ZERO  = SOLID 

C A DATA  CARD  WITH  lOlOlOl  WOULD  INDICATE  4 SOLIDS  AND  THREE  AIR  GAPS 
C IDD  AN  array  of  INTEGERS  THAT  INDICATE  THE  TYPE  OF  SOLIDS  1=  GYPSUM 
C 2=  BRICK  3=  WOOD  SUPPOSE  THERE  WERE  4 SOLIDS  GYPSUM»WOOD» BRICK 
C and  GYPSUM  THEN  DATA  CARD  WOULD  BE  PUNCHED  AS  1321 
C AL  = SOLID  THICKNESS  IN  FEET 
C alpha  = HEAT  DIFU5SI0N  COE^FICIEniT 
C RHO  = DENSITY  OF  SOLID  IN  CUBIT  FEET 

C C=SPECIFIC  HEAT  CAPACITY  0^  SOLID 

C AI  = AIR  GAP  DISTANCE  IN  FEET 

C GA  = heat  release  °ER  UNIT  VOLUME  IN  AIR  GAP 

C OUTPUTS 

C T(IfU)  PRIME  = temperature  AT  T I ^’E ( THETA+DTHfTA ) OF  JTH  MODE  IN 
C ITH  solid  layer 

C TA(I)  = TEMPERATURE  OF  MiD-POiNT  OF  AIR  SPACING  BETWEEN  ITH  AND 
C I+l  LAYER  OF  SOLID  THET^  = TIME 

IMPLICIT  DOUBLE  PRECISI0N( A-Hr 0-W) 

DOUBLE  PRECISION  T1 1 r T22 r T33 » T44 r T66 » T88  »T55 

REAL  TOr  AL»ALPHArRHOrC» AlfGAr  X»Y  rGMl»GM2rGM3 

UIMENSlOiN  T(20r  10)  » AL(20)  » AI  (20)  »TA(20)  »DELTAX(20)  r AM  ( 20  ) » H ( 20  ) » 

1 BH(20) »AK(20»10) r3(20»l0) rFPSLOM(20  ) » RHO ( 20 ) ► C ( 20 ) » 5 ( 20 ) r 

2 P(20)»Rl(20)r  Q(20)fID(30)  » G/' ( 20  ) r ALPHA  ( 20  ) 

DIMENSION  X(500»5) »Y(500»5) »NPMX(5)  » A(2n)»B(2n)r 

1 RHOA(IO) »D(20»10) » IDD(20) f ThETAI ( 20 » 10 ) » ThFTA2 ( 20 » 10 ) 

NA=0 

C THETA  = TIME  IN  HOURS 
theta=o,odo 

C AND  EPN  = EMIS5IYITY 

EP=. 835D0 
EPN=.9D0 

C HH=CHARACTERISTIC  HEAT  TRANSFER  COEFFICIENT 
HH=.15D0 
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C CPA=  SPECIFIC  HEAT  CAPACITY  FOR  AIR 
CPA=0.25D0 

C AMO=  VOLUME  AIR  FLOW  THROUGH  THE  LAYERS^  CUBIC  FEET  PER  HOUR 
AMO=  0,DO 

C ROAl=  DENSITY  OF  nR  IN  L3  P£R  C'JpIC  pOOT 
ROA1=,05DO 

C SIGMA  = aOLTZMAN  STEFAN  CONSTANT 
SIGMA=0.i7i4D“8 

C AK  = HEAT  CONDUCTION  COEFFICIENT 

C READ  IN  DATA  THAT  DESCRIBE  PANELS  AND  READ  IN  OTHER  DATA 

READ  90»NN»N»m^ ITEST»T0»  GMlrGMR»GM3 

kk=nn+m 

READ  91f  (10(1)  rI  = l»KK)* 

READ  91» ( IDD( I) » I=1»NN) 

READ  92» (AL( I) » ALPHA! I)  ? phO ( I ) » C ( I ) » 1=1 » NN) 

IF(M.EQ.O)GO  TO  307 

READ  97»  (Aid)  pGA(I)  »I  = 1»M) 

C DTHETA  = TIME  INCREMENT 

307  DTHETA=1.D0/  720. DO 
DO  301  U=1»N 

DO  301  I=1»NN 

301  T(dU)=T(j 

IF(M.EQ.O)  GO  TO  308 
DO  302  1=1 »M 

302  TA(I)=TO 

308  DO  305  I=1»NN 
H(  I)=HH 

305  EP5L0N(I)=EP 

EP5L0N(NN)=EPN 
EPSL0N(1)=EPN 
DO  104  I=lrNN 

104  DELTAX(I)=AL( I) /FLOAT(N-l) 

DO  105  1=1 »NN 

AM( I)=DELTAX(T )**2/( ALPHA(I)*rTHETA) 
BH(I)=(2.+H(I)^DTHETA)/(RH0(I)*C(I)*DELTAX(I) ) 

105  R1(I)=(2.*DTHETA*SIGMA*EP5L0N(I  ) ) / ( DELTAX ( I ) *RH0 ( I ) *C ( I ) ) 

309  PRINT  100»NN»N» Mr ITEST»T0»  6MlrGM2rGM3 

PRINT  lOlr ( ID( I) r I=lrKK) 

print  lOlr (IDO(I) rI=l,NN) 

1F(M.EQ.Q)G0  TO  310 

PRINT  102r  (Aid)  rGA(I)  rI  = lrM) 

310  PRINT  102r  (ALd)  r ALPHAd)  r RHO  ( I ) » C ( I)  r 1 = 1 r MN  ) 

PRINT  311  r (DELTAXd)  r AM  (1 ) r BH  ( T ) r R1  ( I ) r H ( T ) r 1 = 1 r NN) 
print  102r  (EPSLON(I) r I=lrNN) 

G1  = G,M1 

G2=GM2 

G3=GM3 

T5=T0 

TF=T0 

c Plot  values 

X(lrl)=THETA»?i.n.DO 

Xdr2)=THETA*60.D0 

X(lr3)=THETA*60.D0 

X(lr4)=THETA*60.D0 

X(lr5)=THETA*60.D0 

Y(lrl)=(TF) 

Ydr2)  = (TA(l)  ) 

Ydr3)  = (T(2rN)  ) 
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Y(l»4)=(  TA(2)) 

r { 1 » 5 ) s { T ( i+  r N ) ; 

11=2 

PRINT  37 
PRINT  88 

PRINT  89rTFf  T(  1 »1)  » Ta  ( 1 ) ? T ( 2 ^ k)  ) » TA  ( 2 ) r T ( 3 r 1 ) » TA  ( 3 ) » T ( 4. , N ) ? THETA 
MM=N~1 
1 1 = 1 

IA  = 1 

DO  106  JJ=lnN 
DO  106  LL=1»NM 
106  G(LL» JJ)=O.ODO 

K = 0 

C COMPUTE  TF  FORMULA  7 

IF(THETA.GE.1.D0)G0  TO  415 
TF=TO+101520.DO*THETA/(60.DO*THETA+4.DO) 

GO  TO  417 

415  IF( theta. GE.1.9DO)GO  TO  416 

TF=( 926.00-^42. DO ♦THETA-n.0l31D0*  ( 120 . DO-60 . DO^THETA ) *^2 ) *1 . 3n0 
1 -»-32.D0 

GO  TO  417 

416  TF=(926.D0-»-42.D0*THETA)  ♦1.8D0-J-32.D0 
C COMPUTE  T(l»l)  AND  TPPIME 

1+17  IJ=I0D(1) 

C TEST  PROPERTIES  FOR  CURRENT  SOLID 

CALL  PROP(T(  1 rl)rA<(  1 » 1 ) » THET A 1 ( 1 f 1 ) r THET A2 ( I » 1 ) » THFT A » 

1 DTHETArG(l»l) ,GlrG2»63» IJ) 

924  AKK=AK(1»1) 

TPRIME  =T(1»1)■^R1(1)  ♦(TF-T(1,1)  ) ♦ ( TF  + T ( 1 » 1 ) -^920  . DO  ) ♦ ( (TF-^46^.D0) 

1 ♦♦2-MTdf  l)■^■460.DO)♦♦2)-^BH(l)♦oflBS(TF-T(l,l)  ) ♦♦  . 25*  ( TF-T  ( 1 » 1 ) ) 

2 (G(1»1)*DELTAX(D**2)/(AKK*AM(1)  ) -2.D0/AM(  1)*(T(  1,  1)-T(1»2)  ) 
T(1»1)=TPRIME 

C COMPUTE  TPRIiME(Ifj)  FORVULA  1 
5 IF(AM(I) .LE.2.D0)S0  TO  77 

I J=IOD( I ) 

DO  10  J=2»MM 

C TEST  PROPERTIES  FOR  CURRENT  SOLID 

CALL  PROP(T(  I »J)»AK(  I»J  ) » THETAl ( I » J) » THETA2 ( T » J) » THETA » 

1 DTHETA»G(I» J) »G1»G2»G3»IJ) 

522  AKK=AK(I»J) 

T(I»J)=T11(AM(I) »T(I» J-1) »T(I, J-H) »T(I»J)»G(I»J) rAKK  r 
1 DELTAX(D) 

10  CONTINUE 

IF( I .EQ.NN) GO  TO  45 
IF(  ID(  IA-+I)  .EO.O)  GO  TO  12 
C SOLID  TO  SOLID  COMPUTATION 
IJ=IDD(I ) 

C TEST  PROPERTIES  FOR  CURRENT  SOLID 

CALL  PROP(T(  I »N)»AK(  I»N  ) r THETAl ( I » N ) » THETA2 ( T » N) f THETA » 

1 DTHETA»G(Ir J) rGl»G2»G3r IJ) 

932  AKK=AK(I»N) 

DO  875  JK=1»NN 

A( JK)=(RHO( JK)*C( JK)*dELTAX( JK) )/2.D0 
3( JK)=(DELTAX( JK)*DTHETA)/2,D0 
D( JK»N)=(AK( JK»N)*DTHETA)/DELTAX( JK) 

875  D( JK» 1)=( AK( JK» 1) *DTHETA ) /DELTAX ( JK ) 

T(  I rN)  =T55(  A (I)»T(I»N)»A(I-*-l)»T(I  + l»l)r0Cl)rG(I»N)rB(I-fl)» 

1 G(I  + id)  »D(IfN)  »T(I»N-1)  »T(I-»-l»2)  ) 
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T(I+lf 1)=T{I?N) 

1 = 1+1 

GO  TO  5 

C COMPUTE  TPRIME  FORMULA  2 AIR  SPACE  COMPUTATIOf^ 

12  IF(AM( I ) .LE.2.D0) GO  TO  75 

K=K+1 
IJ=IDO( I) 

C TEST  PROPERTIES  ^OR  CURRENT  SOLID 

CALL  PROP(T(  I »N>»AK(  I»N  ) , THETAl ( I » N ) r THETA2 ( I » N ) » THETA » 

1 OTHETAf  GdoN)  ,G1»G2»G3»IJ) 
bll  AKK=AK(I»N) 

T(IfN)=T22(AiM(I)  fT(Ir\j)».T(IrN-l)»Rl(I)»T(I  + l»l)»BH(I)»TA(K)» 

1 G(I?N) »AKK»DELTAX(I) ) 

C COMPUTE  TAPRImE  FORMULA  3 

RHOA (K) =39. 67400/ (TA ( K) +460.0n) 

P(K)=(H(k)*DTHETA)/(AI (<) *CPA*RH0A(K) ) 

GA(K)-0,000 

b(K)=( AM0^DTHETA)/(AI (K) ) 

U(<)=(H(K+1)*DTHETA)/(AI (K)*CPA+RH0A(K) ) 

GA A=GA ( K ) 

TA(K)=T33(TA(K) >P(K) »T(I»N) » 9 ( K ) » T ( I +1 » 1 ) f S ( K ) » G AA  ) 

C COviPUTE  TPRIME(I  + 1»1)  FOrMJLA  4 
IF(AM( I+i) .LE.2.DO)GO  TO  80 
IJ=IDD(I+1) 

T(I+1p1)=T44(T(I+1»1) fOl(I)»T(T?N)»AM(I+l),RH(I+l)»TA(K)rT(I+l»2)» 
1 G(I  + ld)  »AKKOELTAX(l  + l)  ) 

IF(K.E9.vi)k  = 0 

39  IF(I.EQ.NN)GO  TO  45 

1 = 1 + 1 
1A=IA+1 
GO  TO  5 

C DECIDE  BETWEEN  FORMULAS  5 AND  8 

45  IF( ITEST.NE.1) GO  TO  37 

C COMPUTE  T(NN»N)  FORMULA  6 

IF(AM(NN) .LE.2.D0)30  TO  73 
T1  = TU 

IU=IDj(NN) 

C TEST  PROPERTIES  FOR  CURRENT  SOLID 

Call  prop(T(  nn»n  )»ak(  nn»n) »thetai ( nn»n) »theta2(  nn,n)»theta» 

1 DTHETArG(  NNmN)  »Gl»G2»G3dJ) 

431  AKK=AK(NN»N) 

T(NN»N)=To6(T(NN»N) »AM(NN) »T(NN»N-1) »R1 (NN) »T1 »3H(NN) » 

1 G(NN»N) r AKKf DELTAX(NN) ) 

GO  TO  209 

C compute  T(NNrN)  FORMULA  3 
37  IF(AM(NN) .LE.2.DO)GO  TO  79 

IJ=IQD(NN> 

C TEST  PROPERTIES  POR  CURRENT  SOLID 

Call  prop(t(  nn»n  )»a<(  nn»n) »thetak  nn»n) »theta?{  nn»n) »theta» 

1 0THETA»G(  NN»N)»Gl»G2»G3fIJ) 

432  A<K=AK(NN»N) 

T(NN»N)=T88(T(NN»N)  »T(NN»N-D  » TS » AM ( NN) » BH ( NN) » 

1 G(nN»N) » AKK»DELTAX(NN) »R1(Nn) ) 

P5=0.15DO+DTHETA/(0.5DO*CPA^ROA1) 

TS=TS+PS+DABS(T(NN»N)-T5)*+.25*(T(MNrM)»TS) 

209  IF(M0D(NA» 10) .NE.O)GO  TO  57 

print  89»TFd(l»D»  TA  ( 1 ) » T ( 2 r N)  » TA  (2 ) » T ( 3 » 1 ) » TA  ( 3 ) » T ( 4 » N ) , THETA 
C PLOT  VALUES 
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X(II?1)=THETA»60.D0 
X(  II»2)=THETA«^60.D0 
X(II»3)=THETA*60.DO 
X(II»4)=THETA*60.D0 
X(II»5)=THETA*60.D0 
Ydlf  1)  = (TF) 

Y(II»2)=(TA(1) ) 

Y(II»3)  = (T(2»iM)) 

Y(II»4)=(  TA(2)) 

Y(IIf5)=(T(4»N) ) 

11=11+1 

b7  IFdHETA.GT.  .B33D0)  9O  TO  200 

NA=iMA  + l 

C IMCREMENT  TIME 

theta=theta+jtheta 

GO  TO  1 

200  11=11-1 

PRINT  87 
NROW=500 
NRV|X(1)=II 
NRMX(2) =I I 
NRMX(3)=II 
NRMX(4) =II 
NRMX(5)=II 
NARbS=5 

CALL  PLOTS  (NAogS»  X»  YfNRMXrNROvV) 
STOP 

78  PRINT  81»AM(NN) 

STOP 

79  print  32rAM(NN) 

STOP 

77  PRINT  83»AM(I) 

STOP 

BO  PRINT  84»  AM(T+1) 

STOP 


75  PRINT  85»  AM(I) 

STOP 


81 

F0RMAT( 1X» 

♦ V| 

NOT 

GREATER 

than  2 

FORMULA 

6 

M=*  »E14,8) 

82 

FORMAT (1X» 

♦ M 

not 

greater 

than  2 

FORMULA 

8 

M=»  f E14.8) 

83 

FORMATdXr 

• M 

MOT 

GREATER 

than  2 

FORMULA 

1 

M=’ »E14.8) 

64 

FORMATdXf 

’M 

not 

greater 

than  2 

FORMULA 

4 

M=»  »E14.8) 

65 

FORMATdXf 

f M 

MOT 

greater 

than  2 

FORMULA 

2 

M=*  f E14.8) 

68 

FORMAT ( ’ 

TF 

T(l,l) 

Ta(1) 

1 TA(2) 

T(3» 

1) 

TA(3) 

T(4 

» 

N)  THETA 

101  FORMAT(20I3) 

100  F0RMAT(1X»4I6»5F11.3) 

98  FORMAT(10F6.0) 

92  F0RMAT(  4F6.0) 

97  FORMAT(2Fb.O) 

91  FOR,MAT(20I1) 

216  F0RMAT(lX»lX»2l3»3Fm.2) 

90  FORMAT(4I2»5F6.0) 

89  F0RMAT(lX»Fl2.3»lX»Fl2.3»lXrFl2.3»lX»F12.3»lX»Fl2.3»lX»Fl2.3» 

1 F12.3» 1X»F12.3» 1X»F6.3) 

87  FORMATdHl) 

102  F0RMAT(lXr4Fl4.3) 

312  F0RMATdX»4D26.16) 

311  F0RMAT(1X»F14.8»4026.16) 
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66  FORMAT ( IX » FX2, 3 » IX » Fit* 5 j IX » F12. 3» IX » F12 . 3 ? 13X » F12 , 3 ) 

317  F0RMAT(iX»3F12.3fD2^.I&?3Fl2,3) 

END 


C 

2 


DOUBLE 
IMPLICI 
FORMULA  1 
T11=(T1 
FORMAT ( 


PRECISION  function  Til ( AM » T1 » T2 » T3 » G r AK » DELTAX ) 

T DOUBLE  PRECISI0N(A-H»0-W) 

•»-T2)/AM+  (1«D0-2.D0/AM)*T3  + ( G*DELTAX**2  ) / ( AK*Am 

♦ formula  1») 


) 


RETURN 

END 


c Formula  2 

DOUBLE  PRECISION  FUNCTION  T22 ( AM » T1 r T2 r Rl » T3 » H » TA » G» AX , QELTAX ) 
IMPLICIT  DOUBLE  PREC ISI ON ( A»H » 0-W ) 

T22=Tl+2.D0/Avt*(T2-Tl)”Rl*(Tl-T3)*(Tl+T3+920.D0) ♦( ( Tl+460 . DO ) **2+ 
1 (T3+460.DU) +*2) ”H*DABS(T1“TA)**,'>5*(T1-TA)+(G*DELTAX**2) /( AK^Av) 

1 F0RMAT(lX»3F12,3r026el6p3Fl2.3) 

2 FORiMAT(  » FORMULA  2*  ) 

RETURN 

End 


C FORMULA  3 

DOUBLE  PRECISION  FUNCTION  T33 { TA » P » T1 r 0 » T2 » S r GA ) 

IMPLICIT  DOUBLE  PREC ISI ON ( A-H » 0“W ) 

T33=TA  +P*DAB5(T1“TA)**,25^(T1-TA)-Q^DABS(TA-T2)**.25*(TA-T2)+ 
i S+{T1“T2)+GA 

2 F0RMAT(*  formula  3») 

RETURN 

end 


C 


2 


Formula  4 


1 


function  T44 
IMPLICIT  DOUBLE  PRECISI0N(A-»H 

T44=T1  +R1*(T2»T1)*(T2+T1+920 
”2.D0/AM*(T1-T3) +H*DABS ( TA-T1 
formats  formula  4»r 

RETURN 

End 


T1 » R1 » T2  » AM» 

0“i^) 

DO)*( (T2+460 
**.25*(TA»T1 


HfTA»T3»G»AKoELTAX) 

•D0)*^2+(T1+460.D0  ) ) 

)+(G*DELTAX**2)/(AK*AM) 


C FORMULA  6 

DOUBLE  PRECISION  FUNCTI 0NT66 ( Tl r AM » T2 » R1 » TO r H » G» AK » DELTAX > 
IMPLICIT  DOUBLE  PRECISION! A-H» 0-W) 

T66=Tl+2.D0/AM*(T2-Tl)-Rlt(Tl-T0)^(Tl+T0-^920.D0)*( (T1+460.D0)**2 
1 ^(T0+460.D0)**2)-H*DABS(T1“T0)**.25*(TI-Tn)+(G*DELTAX**2)/(AK*AM) 
2 FORMAT! ♦ FORMULA  6») 

RETURN 

End 
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C FORMULA  8 

DOUBLE  PRECISION!  FUNCTION  Tga  ( T1 » T2 » TS r AM r H f Gr  AK  » OELT AX ? Rl ) 
IMPLICIT  double  PRECI5I0N(A-H»0-W) 

T8a=Tl+2.D0/AM*(T2-Tl)-H*DABS(Tl-TS)**.25  ♦ ( Tl-TS) + (G*DELTAX+*2 ) / 
1 ( AK+AM) 

2 FORMAT(»  FORMULA  8») 

RETURN 

LND 


For  mul  a 

double  PRECISION  FUNCTION  T55  ( Al » T1 » A2  » T2 » B1  r G1 » B2 » G2 » D1 » T3.»  OR » 
1 T4) 

IMPLICIT  DOUBLE  PRECISION  (A-H»0-W) 

T55=  ( A l*Tl  + A2^T2  + Bl’*‘Gl  + B2*G2+Dl*  ( T3-T1 ) -D2*(T2-T4))/(Al  + A?) 

RETURN 

END 


subroutine  PR0P(T1» AK»THETA1,THETA2 » THETA »DDTHET»G»Gm1 6^3 » IU 

1 ) 

IMPLICIT  DOUBLE  PREC I 5l ON ( A-H » 0-W ) 

C TEST  properties  FOR  GYRSUv 
GO  TO  (501»5U2»503) » lU 

501  IF(  Tl.GE.O.ODO.ANO.  T1 . LE . 200 . DO ) AK= . 5D0 

IF(  T1.GE.200.D0.ANO.  Tl . LE . 400 . 00 ) AK= . 3D0 

IF(  T1.GE.400.D0.AND.  T1  . LE. 2000 . DO ) AK  = . 2D0  + T1/4000 . DO 

AK=AK/4.0 

IF(  T1.GE.212.D0.AND.T1.LE.220.D0  ) THETA1=THETA 

IF(  T1.GE.212.D0)THETA2  =THETA-THET A 1 

IF(THETA2.LT.DDTHET)G0  TO  522 

IF(THETA2  . GE . 0 . OOO . AND . THETA2  . LE . 200 . DO+DDTHET ) G =- 

1 (GM1+14100.D0) /(200,dO*DDTHET) 


C 

502 


GO  TO  522 

TEST  properties  BRICK 

IF(  Tl.GE.O.ODO.ANO.  T1 . LE . 200 . DO ) AK=1 . DO 

IF(  T1.GE.200.D0.AND.  T1 . LE. 2000 . DO ) AK  = . 46D0+2 . D-4*T1 

IF{  T1.GE.212.D0.AND.T1.LE.220.D0  ) THETA1=THETA 

IF(  T1.GE.212.D0>THETA2  =THETA-THETA1 

IF(THETA2  . LT . DDTHET ) GO  TO  522 

IF(THETA2  . GE . 0 . OOO . AND . THETA2  .LE.  100 .DO*DDTHET) G =- 


I G;M2/(100.D0*DDTHET) 
GO  TO  522 


C TEST  PROPERTIES  FOR  W00D»D0UGLAS  FIR 

503  AK=.065D0 

IF(  T1.GE.212.D0.AND.T1.LE.220.D0  ) THETA1=THETA 

IF(  T1.3E.212.D0)THETA2  =THETA-THETA1 

IF(THETA2  .LT.DDTHET)G0  to  522 


IF(THETA2  .GE.O,ODO.AND,THETA2  .LE.  200 . DO^ODTHET ) G 

1 GM3/(200.DO*DDTHET) 

522  RETURN 
END 
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